{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|      | 算法 / 技巧                             | 核心思路                                    | 典型用途                   | 代表代码 / 步骤                                      | 误差 / 特性概览                      |\n",
    "| ---- | --------------------------------------- | ------------------------------------------- | -------------------------- | ---------------------------------------------------- | ------------------------------------ |\n",
    "| 1    | 快速逆立方根（Fast Inverse Cubic Root） | 浮点位级“魔数”初值 + 牛顿迭代               | 物理模拟、光照、GPU shader | `i = 0x2a9f7893 - i/3; y *= 4/3 - (x/3)*y²`          | 初值误差约 4%，一次迭代后 <0.1%      |\n",
    "| 2    | Schraudolph 快速 log/exp                | 利用浮点指数偏移，直接线性估算              | 实时机器学习、图形         | `y = i * 2^-23 - 127`（再做多项式校正）              | 粗估 1% 级误差，可经低阶多项式修正   |\n",
    "| 3    | 快速平方根初值（Fast Sqrt Seed）        | 指数右移 + Remez 初值 + 牛顿迭代            | 经典图形流水线、DSP        | `i = (i >> 1) + 0x1fbb4000;` 后 `y = 0.5(y + x/y)`   | 初值误差 <3%，迭代一次逼近 IEEE 精度 |\n",
    "| 4    | De Bruijn 位操作                        | 特殊序列乘法映射到位索引                    | 位运算、压缩、编解码       | `(v & -v)*0x077CB531 >> 27` 查表                     | O(1) 求 ctz/clz，无乘除指令          |\n",
    "| 5    | 整数位扫描开方                          | 按位减法“试探式”构造平方根                  | 固定点、嵌入式 MCU         | 从最高偶数位依次下移，比较减法                       | 无浮点，仅移位加减；复杂度 O(log n)  |\n",
    "| 6    | 快速 atan2 近似                         | 利用象限对称性 + 分段线性                   | FPS、实时图形、导航        | `angle = π/4 - π/4 * r` 等分段表达                   | 误差约 ±0.005 rad，可调节常数控制    |\n",
    "| 7    | Kahan/Neumaier 补偿求和                 | 记录并补偿舍入误差                          | 科学计算、金融统计         | 维护 `c` 补偿量，逐项修正                            | 显著降低长序列累加误差               |\n",
    "| 8    | float-float “伪双精度”                  | 两个 float 打包在 double 尾数               | WebGL、旧 GPU 高精度需求   | 利用 `double` 存 high/low 部分，累计时拆分           | 约可提升 20–30 bit 精度，需小心溢出  |\n",
    "| 9    | SWAR 并行位计数（Hamming Weight）       | 借助掩码与移位，把位分组累加                | 位图处理、压缩、密码学     | 分段掩码 `0x5555...`、`0x3333...` 等逐步求和         | O(1)、无循环，固定指令数             |\n",
    "| 10   | 除常数的魔术乘法（Magic Division）      | 预计算“乘法逆”和位移替代除法                | 编译器优化、数据压缩       | `q = (uint64_t(n) * M) >> s`                         | 精确结果，常数时间，无除法指令       |\n",
    "| 11   | “Magic Bias” 快速 floor/round           | 利用双精度尾数大于 int 精度，添加偏置再减回 | 数值转换、栅格化           | `double y = x + magic; int k = (int)(*(uint64*)&y);` | 适用范围有限，极快，无分支           |\n",
    "| 12   | Bhaskara I 正弦近似                     | 用简单有理式逼近 sin                        | 实时音频、嵌入式图形       | `sin(x) ≈ 16x(π - x)/(5π² - 4x(π - x))`              | 在 [0, π] 误差 ~0.0019，运行时仅乘加 |\n",
    "| 13   | 快速 `pow` 位 hack（Fast Pow Approx）   | 通过 `log2` 的线性近似 + 指数重解释         | 着色器、游戏 AI、特效      | `i = int(b * (i - 1064866805) + 1064866805)`         | 误差可达百分级，胜在速度             |\n",
    "| 14   | Ziggurat 采样算法                       | 分层矩形表 + 快速拒绝采样产生高斯数         | 随机数生成、图形噪声       | 预生成 256 层表；路径判断 + 位操作                   | 每次期望常数时间，统计质量优良       |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这些技巧都带有浓厚的“黑魔法”色彩：充分挖掘浮点 / 整数编码结构，再布上一层迭代或补偿，既快又有趣。即便硬件日益强大，它们仍然是数值算法里的奇妙宝藏。\n",
    "\n",
    "高手潜规则是：\n",
    "\n",
    "1. **充分利用底层表示**（位模式、指数偏移、尾数宽度）。\n",
    "2. **预计算或查表** 把复杂操作前置成本。\n",
    "3. **少量迭代 / 拒绝采样** 取代通用但缓慢的函数调用。\n",
    "\n",
    "掌握它们不仅能在特殊场景下显著提速，也能帮助你更深入理解计算机的数值体系和硬件特性。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 快速逆立方根（Fast Inverse Cubic Root）\n",
    "\n",
    "- **灵感**：延续 Quake III `Q_rsqrt` 的思路，把操作对象换成立方根。\n",
    "- **做法**：将输入浮点数重解释为整数，施加经验“魔数”`0x2a9f7893` 调整指数后，再转回浮点；最后执行一次牛顿迭代。\n",
    "- **用途**：实时光照、物理模拟中需要 `1/∛x` 的场景。\n",
    "- **误差表现**：初始近似约 4%，一次迭代后通常 <0.1%。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "float fast_inv_cbrt(float x) {\n",
    "    float x_third = x / 3.0f;\n",
    "    uint32_t i = *(uint32_t*)&x;\n",
    "    i = 0x2a9f7893 - i / 3;\n",
    "    float y = *(float*)&i;\n",
    "    y = y * (4.0f / 3.0f - x_third * y * y);\n",
    "    return y;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Schraudolph 快速 `log(x)` / `exp(x)`\n",
    "\n",
    "- **核心**：直接利用 IEEE 754 浮点指数偏移量，构建对数的线性近似。\n",
    "- **实现**：将浮点的 32 位模式当整数，乘以 `2^-23` 并减去偏移 127，得到 `log2(x)` 粗估；配合低阶多项式校正尾数提高精度。\n",
    "- **应用**：实时机器学习、渲染曲线拟合、shader 中的快速 log/exp。\n",
    "- **误差**：基础版本约 1%，通过多项式校正可进一步降低。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "float fast_log2(float x) {\n",
    "    union { float f; uint32_t i; } vx = { x };\n",
    "    float y = (float)vx.i * 1.1920928955078125e-7f; // 1 / 2^23\n",
    "    return y - 127.0f;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. 快速平方根初值（Fast Sqrt Seed）\n",
    "\n",
    "- **原理**：平方根相当于指数减半，先对指数右移，再对尾数用 Remez 拟合的常数调整。\n",
    "- **流程**：位级初值后做一次标准牛顿迭代 `y = 0.5 (y + x/y)`，即可达到单精度的误差水平。\n",
    "- **场景**：无硬件 `sqrt` 指令的老平台、音频 DSP、图形流水线。\n",
    "- **特点**：一次迭代就非常精确，成本低。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "float fast_sqrt(float x) {\n",
    "    uint32_t i = *(uint32_t*)&x;\n",
    "    i = (i >> 1) + 0x1fbb4000;\n",
    "    float y = *(float*)&i;\n",
    "    y = 0.5f * (y + x / y);\n",
    "    return y;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. De Bruijn 位操作\n",
    "\n",
    "- **背景**：需要在常数时间内找到最低或最高 set bit 的位置。\n",
    "- **诀窍**：De Bruijn 序列保证 `n` 位子串唯一，利用 `v & -v` 提取最低位，再乘以特定常数并右移，即可索引到正确位置。\n",
    "- **用途**：位图、压缩算法、空间索引结构。\n",
    "- **优点**：纯整数运算，避免循环或分支。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "static const int MultiplyDeBruijnBitPosition[32] = {\n",
    "  0, 1, 28, 2, 29, 14, 24, 3,\n",
    "  30, 22, 20, 15, 25, 17, 4, 8,\n",
    "  31, 27, 13, 23, 21, 19, 16, 7,\n",
    "  26, 12, 18, 6, 11, 5, 10, 9\n",
    "};\n",
    "\n",
    "int ctz32(uint32_t v) {\n",
    "    return MultiplyDeBruijnBitPosition[((v & -v) * 0x077CB531U) >> 27];\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. 整数位扫描开方\n",
    "\n",
    "- **方法**：模仿手算平方根，按位“试探”。从最高偶数位开始，如果当前余数够大就减去 `(res + bit)` 并更新结果，否则仅移位。\n",
    "- **用途**：纯整数环境（MCU、FPGA、固定点 DSP）。\n",
    "- **复杂度**：O(log n)，每轮处理两位二进制位。\n",
    "- **特点**：无需乘除法，极其适合硬件受限场景。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "uint32_t isqrt(uint32_t x) {\n",
    "    uint32_t res = 0;\n",
    "    uint32_t bit = 1u << 30;\n",
    "    while (bit > x) bit >>= 2;\n",
    "    while (bit != 0) {\n",
    "        if (x >= res + bit) {\n",
    "            x -= res + bit;\n",
    "            res = (res >> 1) + bit;\n",
    "        } else {\n",
    "            res >>= 1;\n",
    "        }\n",
    "        bit >>= 2;\n",
    "    }\n",
    "    return res;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. 快速 `atan2` 近似\n",
    "\n",
    "- **思路**：利用 `atan2` 对称性，将平面分象限，再对比例 `r` 用线性组合逼近。\n",
    "- **实现**：常见公式 `angle = π/4 ± π/4 * r`，根据象限调整；对常数微调可控制误差。\n",
    "- **应用**：实时图形、FPS 朝向、导航系统。\n",
    "- **精度**：误差约 ±0.005 弧度，可在速度与精度之间权衡。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "float fast_atan2(float y, float x) {\n",
    "    float abs_y = fabsf(y) + 1e-10f;\n",
    "    float r, angle;\n",
    "    if (x >= 0.0f) {\n",
    "        r = (x - abs_y) / (x + abs_y);\n",
    "        angle = 0.78539816339f - 0.78539816339f * r;\n",
    "    } else {\n",
    "        r = (x + abs_y) / (abs_y - x);\n",
    "        angle = 2.35619449019f - 0.78539816339f * r;\n",
    "    }\n",
    "    return (y < 0.0f) ? -angle : angle;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7. Kahan / Neumaier 补偿求和\n",
    "\n",
    "- **目标**：减轻长序列浮点加法的累积误差。\n",
    "- **机制**：维护补偿变量 `c`，把每次运算丢失的低位再加回来。\n",
    "- **用途**：科学计算、金融分析、统计聚合。\n",
    "- **优势**：实现简单，能显著提高数值稳定性。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "double kahan_sum(const double* a, size_t n) {\n",
    "    double sum = 0.0;\n",
    "    double c = 0.0;\n",
    "    for (size_t i = 0; i < n; ++i) {\n",
    "        double y = a[i] - c;\n",
    "        double t = sum + y;\n",
    "        c = (t - sum) - y;\n",
    "        sum = t;\n",
    "    }\n",
    "    return sum;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "8. float-float “伪双精度”\n",
    "\n",
    "- **背景**：在不支持双精度的环境（如旧 GPU、WebGL）需要更高精度。\n",
    "- **做法**：用两个 `float`（高位 + 低位）表示一个值，通过 Dekker/Two-Sum 等算法确保加法、乘法时把误差分配到低位。\n",
    "- **应用**：高精度地理渲染、科学可视化、复杂矩阵运算。\n",
    "- **优缺点**：可提升有效位 20–30 bit，但实现复杂，需要严格的误差管理。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "9. SWAR 并行位计数（Hamming Weight）\n",
    "\n",
    "- **概念**：SIMD Within A Register（SWAR）技巧，把整数分成多个固定宽度段，同步统计 `1` 的数量。\n",
    "- **操作顺序**：\n",
    "  1. `x -= (x >> 1) & 0x55555555;`\n",
    "  2. `x = (x & 0x33333333) + ((x >> 2) & 0x33333333);`\n",
    "  3. `x = (x + (x >> 4)) & 0x0F0F0F0F;`\n",
    "  4. `x = x * 0x01010101;`\n",
    "  5. `return x >> 24;`\n",
    "- **应用**：快速统计位图人口、CRC、密码学、图像处理。\n",
    "- **特点**：常数指令、无循环或查表，对任何 32 位值都适用，硬件友好。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "10. 除常数的魔术乘法（Magic Division）\n",
    "\n",
    "- **背景**：`n / d`（d 为常数）可转换为 `n * M >> s`，其中 `M` 和 `s` 为预计算常数。\n",
    "- **实现**：\n",
    "  1. 选择 `M ≈ floor(2^(p + s) / d)`（`p` 为字长）。\n",
    "  2. 用 64 位乘法取高位 `>> s`，即得商。\n",
    "- **示例**：`n / 10` 的 `M = 0xCCCCCCCD`, `s = 35`。\n",
    "- **优势**：避免慢速除法指令；现代编译器也使用此策略。\n",
    "- **适用**：只要除数固定即可，无舍入误差（取整结果完全准确）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "11. “Magic Bias” 快速 floor/round\n",
    "\n",
    "- **思路**：双精度 `double` 尾数 52 位，大于 32 位整数。将 `float`/`double` 与一个巨大常数相加再减回，可触发硬件舍入。\n",
    "- **示例**（32 位整型范围）：\n",
    "- **用途**：图形像素定位、音频分桶。\n",
    "- **注意**：仅适用于 |x| < 2^31；需关闭严格别名优化或用 memcpy。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "int fast_floor(float x) {\n",
    "    const double magic = 6755399441055744.0; // 1.5 * 2^52\n",
    "    double y = (double)x + magic;\n",
    "    return ((uint64_t&)y) - 0x8000000000000000ULL;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "12. Bhaskara I 正弦近似\n",
    "\n",
    "- **公式**：对 `x ∈ [0, π]`, $\\sin(x) \\approx \\frac{16x(\\pi - x)}{5\\pi^2 - 4x(\\pi - x)}$\n",
    "- **做法**：先把角度映射到 `[0, π]`，再使用公式；四个乘法、两个加法、一次除法即可。\n",
    "- **误差**：最大绝对误差约 0.0019，对实时音频合成、嵌入式图形很实用。\n",
    "- **延伸**：改良后的 `sin(x)` 公式可与泰勒级数混合获得更小误差。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "13. 快速 `pow` 位 hack（Fast Pow Approx）\n",
    "\n",
    "- **灵感**：沿用 Schraudolph 的 `log2` 思路，把浮点数指数线性化。\n",
    "- **基本近似**：\n",
    "- **含义**：`1064866805` 为 `log2(e)` 相关常数，经线性调整实现 `a^b`。\n",
    "- **优缺点**：速度极快、无调用开销；误差可能高达几%，适用于实时着色器、游戏 AI 或训练中不敏感的启发式。\n",
    "- **改进**：加一次牛顿或多项式校正可显著提升精度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "cpp"
    }
   },
   "outputs": [],
   "source": [
    "float fast_pow(float a, float b) {\n",
    "    int* pa = (int*)&a;\n",
    "    *pa = (int)(b * (*pa - 1064866805) + 1064866805);\n",
    "    return a;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "14. Ziggurat 随机数算法\n",
    "\n",
    "- **目标**：高效生成服从正态（或其他对称分布）的随机数。\n",
    "- **原理**：预构建若干矩形“层”覆盖目标分布（除尾部），每次随机选择一层并在矩形内取样；若落在分布外侧则执行一次简短拒绝采样。\n",
    "- **步骤**：\n",
    "  1. 预计算 256 层的边界 `x[i]` 和 `f(x[i])`。\n",
    "  2. 运行时从均匀随机数中选择层索引和水平坐标。\n",
    "  3. 若点在分布曲线下，直接接受；否则进入尾部处理。\n",
    "- **特性**：期望常数时间，平均只需约 1.02 次循环；适合高性能 RNG、噪声生成、蒙特卡洛模拟。\n",
    "- **魔法感**：通过巧妙几何拆分，大幅减少昂贵开方、对数的调用。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "C++17",
   "language": "C++17",
   "name": "xcpp17"
  },
  "language_info": {
   "name": "C++17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
